/*
 * Main.scala
 *
 * To change this template, choose Tools | Template Manager
 * and open the template in the editor.
 */

package solidutils

import scala.io.Source

object Main {

    /**
     * @param args the command line arguments
     */
    def main(args: Array[String]) :Unit = {
        //        Console.err.println(args.toList.foldLeft("")(_+sep+_))
        args.toList match{
            case "extractFrequencies" ::tail if(tail.size>0)=>
                extractFrequencies(tail(0))
            case "preSNP" :: tail if(tail.size>1)=>
                preSNP(tail(0),tail(1).toInt)
            case _ => Console.err.println(getHelp)
        }
    }

    val getHelp="""Usage solidUtils <cmd> <input>
where <cmd> is:
extractFrequencies - reads a sorted gff-file generated with MaToGff.sh --sort and AnnotageChanges.sh --b and extracts the number of covering starting positions (for each nucleotide) and the frequencies for the nucleotides
    <input> = <gff file>
preSNP - reads a sorted and annotated gff-file generated as in extractFrequencies and extract the number of times nucleotides are present in different starting positions
    <input> = <gff file> <readLength>"""


    val sep:String="\t"
    val nucleotides:List[String]=List("A","C","G","T","N","a","c","g","t","n")
    val complement:Map[Char,Char]=Map(
        ('A'->'T'),('C'->'G'),('G'->'C'),('T'->'A'),('a'->'t'),('c'->'g'),('g'->'c'),
          ('t'->'a'),('N'->'N'),('n'->'n'))

    def preSNP(gffFile:String,readLength:Int):Unit ={
        
        var count:List[readCount]=List()
        var curPos= -1
        var curCount:readCount= new readCount(curPos,readLength)
        for (line <- Source.fromFile(gffFile).getLines){
//            print(line)
            if(!line.startsWith("#")){
                val (id, start, stop, sense, seq)=solidGFFtoObject(line)
                if(curPos== -1){
                    //initialize
                    curPos=start
                    curCount= new readCount(curPos,readLength)
                }
                if(curPos==start){
                    //update count. As the file is sorted, we count all reads in the same position. When that is done, we analyse the result
                    curCount.count(seq, sense)
                }else{
                    count= curCount::count
                    println(curPos+sep+curCount.getSenseCoverage+sep+curCount.getAntisenseCoverage)
                    for(i <- curPos to start -1){
                        //handle result
                        val posCount= new positionCount(i,readLength,Map()++count.filter(_.covers(i)).map((reads)=>((reads.getPosition -> reads))))
                        //println(i+posCount.coverage.slice(0,20).foldLeft("")(_+sep+_))
                        3
                    }
                    //initialise next position
                    curPos=start
                    curCount= new readCount(curPos,readLength) 
                    //remove readCounts that are downstream
                    count= count.filter(_.completelyDownstream(curPos)==false)










//
//
//
//                    for(i <- curPos to start-1){
//                        //print results
//                        var drop= count.keySet.toList
//                        var curCount=newCountStructure() //curCount holds the number of times a letter occurs at a specific position, i.e. each index holds the count from reads starting at the position index positions before the specific position
//                        for (j <-count.keys){
//                            if(j<=i){
//                                val shift=j-i
//                                for (nuc <-nucleotides.map(_.charAt(0))){
//                                    for (strand<- List(true,false)){
//                                        val toAdd= (((count(j))(nuc))(strand)) match {
//                                            case a:List[Int] if(a.length>shift)=>
//                                                drop=drop.filter(_!=j)
//                                                a(shift)
//                                            case _ => 0
//                                        }
//                                        try{
//                                            ((curCount(nuc))(strand))(shift)= toAdd
//                                        }catch{
//                                            case ex:ArrayIndexOutOfBoundsException =>
//                                                val tmpA:Array[Int]=(curCount(nuc))(strand)
//                                                val add:Array[Int]= new Array[Int](shift-tmpA.length+1)
//                                                add(add.length-1)=toAdd
//                                                (curCount(nuc))(strand)=tmpA++add
//                                        }
//                                    }
//                                }
//                            }else{
//                                drop=drop.filter(_!=j)
//                            }
//                        }
//                        //remove keys not needed
//                        for (key <- drop){
//                            count=count-key
//                        }
//                        //analyse
//                        //calculate sum for different positions
//                        val sum = getSumOfCountStructure(curCount)
//                        //nullify outliers
//                        //                        for (nuc <- nucleotides.map(_.charAt(0))){
//                        //                            for (strand <- List(true,false)){
//                        //                                val tmpA:Array[Int]= ((curCount(nuc))(strand))
//                        //                                val E= expectation(tmpA)
//                        //                                (curCount(nuc))(strand) = nullifyOutliers(tmpA, E, standardDeviation(tmpA, E))
//                        //                            }
//                        //                        }
//
//                        //print
//                        print(i)
//                        for (nuc <- nucleotides.map(_.charAt(0))){
//                            for (strand <-List(true,false)){
//                                print(sep+nuc+strand)
//                                for (count:Int <- (curCount(nuc))(strand)){
//                                    print(sep+count)
//                                }
//                            }
//                        }
//
//                    }
                    
                }
            }
        }
    }

    """Sets all positions that deviates from the expected value with more than 3.5*StandardDeviation to zero"""
    def nullifyOutliers(data:Array[Int],Expectation:Double,StandardDeviation:Double):Array[Int] = data.map((x)=>if(Math.abs(Expectation-x)<3.5*StandardDeviation) x else 0)

    def expectation(data:Array[Int]):Double=data.foldLeft(0.0)(_+_)/data.length

    def variance(data:Array[Int],Expectation:Double):Double= (data.map(Expectation-_).map(Math.pow(_, 2)).foldLeft(0.0)(_+_))/(data.length-1)

    def variance(data:Array[Int]):Double=variance(data,expectation(data))

    def standardDeviation(data:Array[Int],Expectation:Double)=Math.sqrt(variance(data,Expectation))

    def standardDeviation(data:Array[Int]):Double=standardDeviation(data,expectation(data))

    def solidGFFtoObject(line:String):(String,Int,Int,Boolean,String) ={
        val l=line.slice(0,line.length-1).split("\t")
        if(l.length!=9){
            throw new Exception("bad line: "+line)
        }
        val id=l(0)
        val start=l(3).toInt
        val stop=l(4).toInt
        val sense= l(6).equals("+")
        if(!(l(6).equals("-")||l(6).equals("+"))){
            throw new Exception("Strange line:\n"+line)
        }
        try{
            l(8).split("\t")(0).split("=")(1)
        }catch{
            case ex:Exception => Console.err.println(l(8))
        }

        val seq= if(sense){
            l(8).split(";")(0).split("=")(1)
        }else{
            l(8).split(";")(0).split("=")(1).map(complement.get(_).get).reverse.foldLeft("")(_+_)
        }
        (id,start,stop,sense,seq)
    }

    def extractFrequencies(gffFile:String):Unit ={
        //        val nucleotides:List[String]=List("A","C","G","T","N","a","c","g","t","n")
        //        val complement:Map[Char,Char]=Map(
        //            ('A'->'T'),('C'->'G'),('G'->'C'),('T'->'A'),('a'->'t'),('c'->'g'),('g'->'c'),
        //              ('t'->'a'),('N'->'N'),('n'->'n'))
        var count:Map[Int,Map[String,Int]]=Map()
        var startPoints:Map[Int,Map[String,List[Int]]]=Map()
        var Min:Int=Integer.MAX_VALUE
        var Max:Int=Integer.MIN_VALUE
        var linenr=0
        println(nucleotides.foldLeft(nucleotides.foldLeft("pos"+sep+"unique startPoints")(_+sep+_)+sep+"coverage")(_+sep+_))
        for ( line <- Source.fromFile(gffFile).getLines){
            linenr+=1
            if(linenr%100000==0)Console.err.println(linenr)
            if(!line.startsWith("#")){
                val l=line.slice(0,line.length-1).split("\t")
                if(l.length==9){
                    //                    Console.err.println("ett")
                    val start=l(3).toInt
                    val stop=l(4).toInt
                    Min= Min.min(start)
                    Max= Max.max(stop)
                    val sense:Boolean=l(6).equals("+")
                    //                    Console.err.println("två")
                    if(!(l(6).equals("-")||l(6).equals("+"))){
                        Console.err.println("Strange line:\n"+line)
                    }
                    val seq= if(sense) {
                        l(8).split("\t")(0).split("=")(1) 
                    }else{
                        //reverse complement
                        l(8).split(";")(0).split("=")(1).map(complement.get(_).get).reverse.foldLeft("")(_+_)
                    }
                    for(i<-0 to seq.length-1){
                        val s:Int=if(sense) (i+1) else -(i+1)
                        val current=seq.charAt(i)+""
                        count.get(start+i) match{
                            case Some(old)=> 
                                old.get(current) match{
                                    case Some(c)=> count+= (start+i ->old.update(current,c+1))
                                    case None => count+= (start+i ->old.update(current,1))
                                }
                                def updateStartPoint(key:String):Map[String,List[Int]]= {

                                    val startPoint:Map[String,List[Int]] = startPoints.get(start+i) match{
                                        case Some(map) => map
                                        case _ =>
                                            throw new Exception("blllllääää")
                                    }
                                    startPoint.get(key) match{
                                        case Some(list) if(!(list.contains(s))) =>
                                            //                                            Console.err.println("irk")
                                            startPoint.update(key, s::list)
                                        case Some(list) => startPoint
                                        case None => startPoint.update(key,List(s))
                                    }
                                }
                                startPoints+= (start+i-> updateStartPoint("tot"))

                                startPoints+= (start+i-> updateStartPoint(current))

                            case None =>
                                count+=(start+i ->Map((current->1)))
                                startPoints+= (start+i->Map( (current->List[Int](s)), ("tot"->List[Int](s)) ))
                                
                        }
                    }
                    //flush old data
                    for(i<-Min to start-1){
                        val sp:List[Int]= startPoints.get(i) match{
                            case Some(startPoint)=>
                                startPoint.getOrElse("tot", List()).size :: nucleotides.map(startPoint.getOrElse(_,List()).size)
                            case None =>
                                0::nucleotides.map((_)=>0)
                        }
                        val cov:List[Int]=count.get(i) match{
                            case Some(c) => nucleotides.map(c.getOrElse(_,0))
                            case None => nucleotides.map((_)=>0)
                        }
                        println(i+sp.foldLeft("")(_+sep+_)+sep+cov.foldLeft(0)(_+_)+cov.foldLeft("")(_+sep+_))
                        startPoints-=i
                        count-=i
                    }
                    Min=start
                }
            }
            
        }

        //print
        for (i <- Min to Max){
            val sp:List[Int]= startPoints.get(i) match{
                case Some(startPoint)=>
                    startPoint.getOrElse("tot", List()).size :: nucleotides.map(startPoint.getOrElse(_,List()).size)
                case None =>
                    0::nucleotides.map((_)=>0)
            }
            val cov:List[Int]=count.get(i) match{
                case Some(c) => nucleotides.map(c.getOrElse(_,0))
                case None => nucleotides.map((_)=>0)
            }
            println(i+sp.foldLeft("")(_+sep+_)+sep+cov.foldLeft(0)(_+_)+cov.foldLeft("")(_+sep+_))
        }
    }


}
